//function  least-squares method for calculate temperature value and density value
#include "TomsonScattering.h"

#include <vector>
#include <algorithm>
#include <functional>
#include <TCanvas.h>
#include <TROOT.h>
#include <TStyle.h>

#include <iostream>

#include <cmath>
#include "Exceptions.h"
#include "sqr.h"

void TomsonScattering::Calculate(CalibData &calibrations, bool useStatWeight, bool useFilters, bool timeSync, bool autoChooseChanels)
{
    m_errors.clear();
    CheckConfig();
    std::vector<ValueWithError> integralPeakNoCalib;
    CalculateSignalIntegrals(integralPeakNoCalib, useFilters, timeSync, autoChooseChanels);
    
    if ( integralPeakNoCalib.size() < m_chanelEnabled.size() + 1 )
        throw ProcessingException(std::string( "Not enough signal data in this shot" ));

    CalculateShotParams(calibrations, integralPeakNoCalib, useStatWeight);
    CleanUp();
}

void TomsonScattering::CalculateSignalIntegrals(std::vector<ValueWithError> &integralPeakNoCalib, bool useFilters, bool timeSync, bool autoChooseChanels )
{
    long int maxI = 1L << m_qbitsADT;
    double minY = 0.0;
    double maxY = 0.0;

    double timeQuant = 0.0;
    double amplification = 1.0;

    // used only for output to std::out
    std::vector<double> amplifacations;
    for(int j = 1; j <= m_chanelNum; j++)
    {
        if ( j > m_chanelEnabled.size () && j != m_chanelNum )
            continue;

        if ( j <= m_chanelEnabled.size() )
        {
            if ( m_chanelEnabled[j-1] == false )
            {
                integralPeakNoCalib.push_back( ValueWithError() );
                continue;
            }
        }

        TH1F * signalHistogram = GetSignalData( j, maxI, timeQuant, minY, maxY, amplification );
        amplifacations.push_back(amplification);

        if( signalHistogram == NULL)
        {
            integralPeakNoCalib.push_back( ValueWithError() );
            continue;
        }

        std::vector<double> signal;

        ReadSignal(signalHistogram, signal);

        double signalLength = signalHistogram->GetNbinsX()*timeQuant;

        if( m_visHelper != NULL )
        {
            if( j <= m_chanelEnabled.size() || j == m_chanelNum )
                m_visHelper->DrawSignal(signal, maxY, minY, 0, signalLength, j, false);
        }

        if ( useFilters )
            MedianLowpassFilter(signal);

        if( m_visHelper != NULL )
        {
            if( j <= m_chanelEnabled.size() || j == m_chanelNum )
            {
                if( useFilters )
                {
                    m_visHelper->DrawSignal(signal, maxY, minY, 0, signalLength, j, true);
                }
                else
                    m_visHelper->DrawEmpty();
            }
        }

        int startInd = 0, stopInd = 0;

        if( j == m_chanelNum )
        {
            startInd = std::max(int(m_signal12Start/timeQuant), 0);
            stopInd = std::min(int(m_signal12End/timeQuant), int(signal.size()) );

            timeSync = false;
        }
        else
        {
            startInd = std::max(int(m_signalStart/timeQuant), 0);
            stopInd = std::min(int(m_signalEnd/timeQuant), int(signal.size()));
        }
        // erase signal data we do not need (second peak too)
        std::vector<double>::iterator it = signal.begin();
        std::advance(it, startInd+1);
        signal.erase(signal.begin(), it);

        if(timeSync)
            stopInd = signal.size();
        else
            stopInd -= startInd;

        int maxIndex = FindMaxValueIndex(signal, 0, stopInd, timeSync);

        if( signal[maxIndex] >= maxY+minY )
            m_errors << "Chanel " << j << ": peak value (" << signal[maxIndex] << ") exceeds maximum (" << maxY+minY << ")" << std::endl;

        if ( autoChooseChanels )
            CheckChanel(signal, maxIndex, j-1);

        ValueWithError backgroundNoise;
        CalcAverageAndError( maxIndex - m_beforeMaxStart, maxIndex - m_beforeMaxEnd + 1, signal, backgroundNoise);

        IntegrateSignal(signal, backgroundNoise, amplification);

        ValueWithError noiseBeforeSignal;
        CalcAverageAndError( maxIndex - m_beforeMaxStart, maxIndex - m_beforeMaxEnd + 1, signal, noiseBeforeSignal);

        ValueWithError noiseAfterSignal;
        CalcAverageAndError( maxIndex + m_afterMaxStart, maxIndex + m_afterMaxEnd + 1, signal, noiseAfterSignal);

        integralPeakNoCalib.push_back( noiseAfterSignal - noiseBeforeSignal );
    }

    if( m_visHelper != NULL )
    {
        m_visHelper->PrintVector(std::string("Integrals without calibration"), integralPeakNoCalib);
        m_visHelper->PrintVector(std::string("Signal amplifications"), amplifacations);
    }

    if ( integralPeakNoCalib.back().val > 0)
    {
        for( int i = 0; i < integralPeakNoCalib.size()-1; ++i )
            integralPeakNoCalib[i].val /= integralPeakNoCalib.back().val;
    }
 };

void TomsonScattering::CalculateShotParams(CalibData &calibData, std::vector<ValueWithError > &integralNoCalib, bool useStatWeight)
 {
    // small
    double r = calibData.smallCalib.radius; // radius
    double l = calibData.smallCalib.distanceAxisLaserRay; //60.200001; // distance
    double Te_coef = calibData.smallCalib.temperatureCoef;
    double teta = calibData.smallCalib.angle;//86;
    double ne_coeff = calibData.smallCalib.densityCoef; //3.84;
    double energy_coef = calibData.smallCalib.energyCoef;

    // background (fon)
    std::vector<ValueWithError> bgCalib(calibData.bgCalib.begin(), calibData.bgCalib.end());

    // spectral
    std::vector<double> chanelCenter(calibData.spectralCalib.specLine.begin(), calibData.spectralCalib.specLine.end());// mid of chanels
    std::vector<ValueWithError> spectralCalib(calibData.spectralCalib.specCoef.begin(), calibData.spectralCalib.specCoef.end());//resultati spectralnoi kalibrovki

    // peak ratio
    std::vector<ValueWithError> peakRatio(calibData.peakRatioCalib.begin(), calibData.peakRatioCalib.end());//resultati otnositelnoi kalibrovki

    std::vector<double> integralPeak; //integration of every chanel with callibrations (of 5).  1ii pik
    std::vector<double> logIntegralPeak;//log of integration of every chanel
    std::vector<double> statWeight;//statvesa

    for(int j = 0; j < m_chanelEnabled.size(); ++j)
    {
        integralPeak.push_back( (integralNoCalib[j].val-bgCalib[j].val)/(peakRatio[j].val*spectralCalib[j].val) );
        logIntegralPeak.push_back( log(integralPeak.back()) );

        double sqrSignal = sqr(integralNoCalib[j].error / ( (integralNoCalib[j].val - bgCalib[j].val) * integralNoCalib.back().val ));
        double sqrBg = sqr( bgCalib[j].error / ( integralNoCalib[j].val - bgCalib[j].val ) );
        double quadSignalBg = sqr( integralNoCalib[j].val * integralNoCalib.back().val * integralNoCalib.back().error / ( (integralNoCalib[j].val - bgCalib[j].val) * sqr(integralNoCalib.back().val) ) );
        double pkRatio = sqr( peakRatio[j].error / peakRatio[j].val );
        double sqrSpec = sqr( spectralCalib[j].error / spectralCalib[j].val );

        statWeight.push_back( fabs(logIntegralPeak.back()) / std::sqrt( sqrSignal + sqrBg + quadSignalBg + pkRatio + sqrSpec) );
    }

    if( m_visHelper != NULL )
    {
        m_visHelper->PrintVector(std::string ("statWeight"), statWeight);

        m_visHelper->PrintVector(std::string ("Logarithms of calibrated integrals"), logIntegralPeak);
        double energy = integralNoCalib.back().val * energy_coef;
        m_visHelper->PrintValue(std::string ("Energy"), energy);

        std::vector<double> invSNR;
        for( int i = 0; i < m_chanelEnabled.size() && i < integralNoCalib.size(); ++i )
            invSNR.push_back( 100.0 * bgCalib[i].val / integralNoCalib[i].val );

        m_visHelper->PrintVector(std::string ("Noise to signal [%]"), invSNR);
    }

    double tetaRad = (3.14159265 * teta) / (360.0) + atan(r/l) / 2.0;//perevod v radian
    double alpha = 1.0/(Te_coef*sqr(sin(tetaRad)));//coefficient for te

    CalcParams(m_Te, m_Ne, m_line, m_points, alpha, chanelCenter, logIntegralPeak, statWeight, ne_coeff, useStatWeight);

    std::vector<double> logIntegralPeakCAL;//log of 1i pik . dlya osi y s shirinoi
    CalcIntegralCorrection(logIntegralPeakCAL, tetaRad, chanelCenter, integralPeak);

    CalcParams(m_Te_chanelWidth, m_Ne_chanelWidth, m_line_chanelWidth, m_points_chanelWidth, alpha, chanelCenter, logIntegralPeakCAL, statWeight, ne_coeff, useStatWeight);
};

void TomsonScattering::CalculateBgCalib(std::vector<ValueWithError>& bgCalib, bool useFilters, bool timeSync)
{
    bgCalib.clear();
    CalculateSignalIntegrals(bgCalib, useFilters, timeSync, false);
    bgCalib.pop_back();
}

void Swap(double &a, double &b)
{
    double t = a;
    a = b;
    b = t;
}

void TomsonScattering::CalculatePeakRatioCalib(std::vector<ValueWithError>& peakCalib, calib_options opts, bool useFilters)
{
    void Swap(double &a, double &b);
    
    if ( opts == peak_normSecond )
    {
        Swap(m_signalStart, m_secondPeakStart);
        Swap(m_signalEnd, m_secondPeakEnd);

        Swap (m_maxOffset, m_secPeakOffset);

        Swap(m_beforeMaxStart, m_beforeSecondPeakStart);
        Swap(m_beforeMaxEnd, m_beforeSecondPeakEnd);
        Swap(m_afterMaxStart, m_afterSecondPeakStart);
        Swap(m_afterMaxEnd, m_aftereSecondPeakEnd);
    }

    peakCalib.clear();
    CalculateSignalIntegrals(peakCalib, useFilters, false, false);
    peakCalib.pop_back();

    if ( opts == peak_normSecond )
    {
        Swap(m_secondPeakStart, m_signalStart);
        Swap(m_secondPeakEnd, m_signalEnd);

        Swap (m_secPeakOffset, m_maxOffset);
        
        Swap(m_beforeSecondPeakStart, m_beforeMaxStart);
        Swap(m_beforeSecondPeakEnd, m_beforeMaxEnd);
        Swap(m_afterSecondPeakStart, m_afterMaxStart);
        Swap(m_aftereSecondPeakEnd, m_afterMaxEnd);
        
        opts = peak_normCustomCh;
    }

    NormPeakRatioCalib(peakCalib, opts);
}

void TomsonScattering::NormPeakRatioCalib(std::vector<ValueWithError>& peakCalib, calib_options opts)
{
    switch(opts)
    {
        case peak_normCustomCh:
        {
            if( peakCalib.size() > m_peakCalibChanel )
            {
                ValueWithError denom = peakCalib[m_peakCalibChanel];
                for( int i = 0; i < peakCalib.size(); ++i )
                    peakCalib[i] = peakCalib[i] / denom;
            }

            break;
        }
        case peak_normSum:
        {
            double sum = 0;
            for( int i = 0; i < peakCalib.size(); ++i )
                sum += peakCalib[i].val;
            for( int i = 0; i < peakCalib.size(); ++i )
                peakCalib[i].val /= sum;

            break;
        }
        case peak_additionalNorm:
        {
            // TODO: not implemented in hardware
        }
    }
}

void TomsonScattering::CheckConfig()
{
    if( m_deltaL.size() != m_chanelEnabled.size() )
    {
        m_errors << "Delta L must be equal to the number of signals" << std::endl;
        throw ProcessingException( std::string("Delta L must be equal to the number of signals") );
    }

    if( m_arc == NULL )
    {
        m_errors << "Supplied archive is not valid" << std::endl;
        throw ProcessingException( std::string("Supplied archive is not valid"));
    }
}

void TomsonScattering::CleanUp()
{
    m_arc = NULL;
    m_visHelper = NULL;
}

void TomsonScattering::LeastSquaresMethod(std::vector<double> &chanelCenter, std::vector<double> &logIntegralPeak, std::vector<double> &statWeight, std::vector<double> &wavelengthDelta, ValueWithError *linePars)
{
    wavelengthDelta.assign( m_chanelEnabled.size(), 0 );
    double xsr = 0, ysr = 0, summstatves = 0, Dw = 0, summdi = 0;
    int i;

    for(i = 0; i < m_chanelEnabled.size(); ++i)
    {
        wavelengthDelta[i]  = sqr( m_laserWavelength - chanelCenter[i] );
        
        if( m_chanelEnabled[i] )
        {
            xsr += wavelengthDelta[i] * statWeight[i];
            ysr += logIntegralPeak[i] * statWeight[i];
            summstatves += statWeight[i];
        }
    }

    xsr /= summstatves;
    ysr /= summstatves;

    for(i = 0 ;  i < m_chanelEnabled.size(); ++i)
    {
        if( m_chanelEnabled[i] )
        {
            linePars[0].val = linePars[0].val + statWeight[i] * (wavelengthDelta[i] - xsr) * logIntegralPeak[i];
            Dw += statWeight[i] * sqr( wavelengthDelta[i] - xsr );
        }
    }

    // y = ax+b
    linePars[0].val /= Dw; // a
    linePars[1].val = ysr - linePars[0].val * xsr; // b

    double d = 0.0;
    for( i = 0; i < m_chanelEnabled.size(); ++i )
    {
        if( m_chanelEnabled[i] )
        {
            d =  logIntegralPeak[i] - wavelengthDelta[i] * linePars[0].val - linePars[1].val;
            summdi += statWeight[i] * sqr(d);
        }
    }

    linePars[0].error = std::sqrt((1.0/Dw * summdi/4.0));//delta a
    linePars[1].error = std::sqrt((1.0/summstatves + sqr(xsr)/Dw)*(summdi/(4.0)));//delta b

    return;
}

void TomsonScattering::BubbleSort(double* tmp)
{
    double trash = 0;
    
    for (int i=1;  i<=3  ;  i++)
    {
        for (int j=0;  j<=1;  j++)     
        {
        
            if (*(tmp + j) > *(tmp + (j +1)))     
            {
                trash = *(tmp + j);
                *(tmp + j) = *(tmp + (j + 1));
                *(tmp + (j + 1)) = trash;
            }
        }
    }
}


// filteres signal with median filter and smoothes signal with low-pass filter
int TomsonScattering::MedianLowpassFilter(std::vector<double> &signal, unsigned int windowSize, double lowPassCoef)
{
    if(signal.size() < windowSize )
        return 0;
    
    double *window = new double[windowSize*sizeof(double)];

    double pred_value = 0.0; // previous smoothed value of signal for low-pass filter
    double coeffi = 0.6;

	std::vector<double> signal_tmp;
	signal_tmp.assign(signal.begin(), signal.end());
	
	for(int j = (windowSize-1)/2; j < signal.size() - (windowSize-1)/2; ++j)
	{
		for(int i = 0; i < windowSize; ++i)
		{
			window[i] = signal_tmp[j+i-(windowSize-1)/2];
		}

		BubbleSort(window);

		pred_value = coeffi*window[(windowSize-1)/2] - (coeffi - 1.0)*pred_value;
		signal[j] = pred_value;
	}

	delete[] window;

	return 0;
}

void TomsonScattering::IntegrateSignal(std::vector<double> &signal, ValueWithError noiseLevel, double amplification)
{
	double sum = 0;
	//integrate the signal in a given interval, minus dust, use amplification

	for( int i = 0; i < signal.size(); ++i)
	{
		sum += (signal[i] - noiseLevel.val) / amplification;
		signal[i] = sum;
	}
}

// returns signal value avaraged in [start, stop) interval and averaging error
void TomsonScattering::CalcAverageAndError(int start, int stop, std::vector<double> &signal, ValueWithError &average)
{
    if ( start < 0 || start >= signal.size() )
    {
        m_errors << "Integrating interval start position out of range. Setting it 0." << std::endl;

        start = 0;
        
        //throw ProcessingException( std::string("Maximim steoffs out of range") );
    }

    if(  stop < 0 || stop >= signal.size() )
    {
        m_errors << "Integrating interval stop position out of range. Setting it " << signal.size()-1 << "." << std::endl;
        stop = signal.size()-1;
    }

    double size = stop - start;
    for (int i = start; i < stop; ++i)
    {
        average.val += signal[i];
    }

    average.val /= size;

    for (int i = start; i < stop; ++i)
    {
        average.error += sqr( average.val - signal[i] );
    }

    average.error = std::sqrt( average.error/size );
}

void TomsonScattering::ReadSignal(TH1F* h, std::vector<double> &signal)
{
    for( int i = 0; i < h->GetNbinsX(); ++i )
    {
	signal.push_back( h->GetBinContent(i) );
    }
}

int TomsonScattering::FindMaxValueIndex(std::vector<double> &signal, int start, int end, bool timeSync)
{
    int maxIdx = start;
    if ( !timeSync )
    {
        for( int i = start+1; i < end; ++i )
        {
            if (signal[i] > signal[maxIdx])
            {
                maxIdx = i;
            }
        }
    }
    else
    {
        if ( m_syncMaxIdx < 0 )
        {
            double timeQuant, minY, maxY, amplification;
            long int maxI = 1L << m_qbitsADT;
            TH1F * signalHistogram = GetSignalData( m_chanelNum, maxI, timeQuant, minY, maxY, amplification );

            std::vector<double> signal12;
            ReadSignal(signalHistogram, signal12);

            int startInd = m_signal12Start/timeQuant;

            // erase signal data we do not need (second peak too)
            std::vector<double>::iterator it = signal12.begin();
            std::advance(it, startInd+1);
            signal12.erase(signal12.begin(), it);

            int maxIndex = FindMaxValueIndex(signal12, 0, signal12.size(), false);
            m_syncMaxIdx = (m_signal12Start+maxIndex*timeQuant + m_maxOffset - m_signalStart)/timeQuant;
        }

        maxIdx = (m_syncMaxIdx >= 0)?m_syncMaxIdx:0;
    }

    return maxIdx;
}

TH1F* TomsonScattering::GetSignalData( int chanelNum, long int maxI, double &tquant, double &minY, double &maxY, double &amplification )
{
	TSignal* signalData = m_arc->Signal(chanelNum);
        amplification = m_arc->Amplification(chanelNum);

	if(signalData == NULL)
            return NULL;
		
	minY = signalData->GetYShift();
        double vquant = signalData->GetYQuant(); // in Volts
        std::cout << " vquant: " << vquant << std::endl;
	maxY = minY + maxI * vquant;
	
	tquant = signalData->GetXQuant(); // in mks

	std::cout << " tquant: " << tquant << std::endl;

	double t0 = signalData->GetXShift(); // in mks
	unsigned int size = signalData->GetSize();

	char hName[100];
	sprintf(hName,"%s_hist",signalData->GetName());
	TH1F* h=(TH1F*)(signalData->CreateHistogram()->Clone(hName));

	delete signalData;

	return h;
}

bool TomsonScattering::CheckChanel(std::vector<double> &signal, int maxIdx, int chanelNum)
{
    int noiseMaxIndex = FindMaxValueIndex(signal, 0, maxIdx - m_beforeMaxEnd, false);
    if ( ( signal[noiseMaxIndex] >= signal[maxIdx]*0.75f ) && ( chanelNum < m_chanelEnabled.size() ) )
    {
        m_chanelEnabled[chanelNum] = false;
        m_errors << "Chanel " << chanelNum + 1 << " was switched off: signal to noise validation" << std::endl;
        return false;
    }

    return true;
}

// the main function - here is the core of the whole program
void TomsonScattering::CalcParams(ValueWithError &Te, ValueWithError &Ne, std::vector<std::pair<double, double> > &line,
        std::vector<std::pair<double, ValueWithError> > &points, double alpha, std::vector<double> &chanelCenter,
        std::vector<double> &logIntegralPeak, std::vector<double> &statWeight, double Ne_coef, bool useStatWeight)
{
    const double e = 2.718281828459045;

    std::vector<double> wavelengthDelta;
    ValueWithError lineParams[2];//parameters of directly is counted by the method of least squares a[0] = a, a[1] = b, a[2] = delta_a, a[3] = delta_b in ax+b = y
    
    std::vector<double> weights;
    if( useStatWeight )
        weights.assign(statWeight.begin(), statWeight.end());
    else
        weights.assign(statWeight.size(), 1);

    LeastSquaresMethod(chanelCenter, logIntegralPeak, weights, wavelengthDelta, lineParams);

    Te.val = -alpha / lineParams[0].val;
    Te.error = alpha * lineParams[0].error / sqr( lineParams[0].val );
    Ne.val = Ne_coef * pow( e, lineParams[1].val ) / std::sqrt( -lineParams[0].val );

    Ne.error = Ne_coef*pow(e, lineParams[1].val) * std::sqrt( sqr(lineParams[1].error) / fabs(lineParams[0].val)
		+ sqr(lineParams[0].error / (2.0 * pow(fabs(lineParams[0].val), 1.5))) );

    line.clear();
    line.push_back( std::make_pair(0, lineParams[1].val) );

    line.push_back( std::make_pair(wavelengthDelta[4], lineParams[0].val*wavelengthDelta[4] + lineParams[1].val) );

    points.clear();
    for( int i = 0; i < m_chanelEnabled.size(); ++i )
    {
        if( m_chanelEnabled[i] )
            points.push_back(std::make_pair(wavelengthDelta[i], ValueWithError(logIntegralPeak[i], fabs(logIntegralPeak[i])/statWeight[i])));
    }
}

void TomsonScattering::CalcIntegralCorrection(std::vector<double> &logIntegralPeakCAL, double tetaRad, std::vector<double> chanelCenter, std::vector<double> integralPeak)
{
    //calculate the temperatute use width of chanels
    double sigmaExp = std::sqrt((m_Te.val*8.7))*sin(tetaRad);//sigma - parametr to construct gauss curve for obtaining temperature

    //integration of gauss part in chanels
    for( int i = 0; i < m_chanelEnabled.size(); ++i )
    {
        double gausSum = 0;//integration of an ideal Gauss in these intervals
        for ( double index = chanelCenter[i] - m_deltaL[i]; index < chanelCenter[i] + m_deltaL[i]; index+=0.005 ) //integration po shirine of chanel
        {
            gausSum += (1.0/(2.0*sigmaExp))*exp(-sqr(index - m_laserWavelength)/(2.0*sqr(sigmaExp)));
        }

        gausSum /= (1.0 + 2.0*m_deltaL[i]/0.005);

        double cal = gausSum / ((1.0/(2.0*sigmaExp))*exp(-sqr(chanelCenter[i]-m_laserWavelength)/(2.0*sqr(sigmaExp))));

        logIntegralPeakCAL.push_back( log(integralPeak[i] / cal) );
    }
}

 

ValueWithError TomsonScattering::GetTe(bool chanelWidth)
{
    if( chanelWidth )
        return m_Te_chanelWidth;
    else
        return m_Te;
}

ValueWithError TomsonScattering::GetNe(bool chanelWidth)
{
    if( chanelWidth )
        return m_Ne_chanelWidth;
    else
        return m_Ne;
}

std::vector<std::pair<double, double> > TomsonScattering::GetLine(bool chanelWidth)
{
    if( chanelWidth )
        return m_line_chanelWidth;
    else
        return m_line;
}

std::vector<std::pair<double, ValueWithError> > TomsonScattering::GetPoints(bool chanelWidth)
{
    if( chanelWidth )
        return m_points_chanelWidth;
    else
        return m_points;
}
